﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using RSPJ.Struct;
using RSPJ.Windows;

namespace RSPJ.Function
{
    public class funcStretch
    {
        // 任意拉伸
        public Files LinearStretch(Files file, int a, int b, int c, int d)
        {
            Files outFile = file.Clone();
            outFile.Name += "_Stretched";

            for (int i = 0; i < file.Head.Bands; i++)
            {
                for (int j = 0; j < file.Head.Lines; j++)
                {
                    for (int k = 0; k < file.Head.Samples; k++)
                    {
                        if (file.IntArray[i, j, k] < a)
                        {
                            outFile.IntArray[i, j, k] = c;
                        }

                        else if (file.IntArray[i, j, k] > b)
                        {
                            outFile.IntArray[i, j, k] = d;
                        }
                        else
                        {
                            outFile.IntArray[i, j, k] = (int)((d - c) * (file.IntArray[i, j, k] - a) / (b - a) + c);
                        }

                    }
                }
            }
            return outFile;
        }

        // 2%线性拉伸
        public Files twoLinearStretch(Files file)
        {
            Files outFile = file.Clone();
            outFile.Name += "_2%Stretched";

            //计算灰度直方图
            int[,] Histo = new int[7, 256];

            for (int i = 0; i < file.Head.Bands; i++)
                for (int j = 0; j < file.Head.Lines; j++)
                    for (int k = 0; k < file.Head.Samples; k++)
                    {
                        for (int p = 0; p < 256; p++)
                        {
                            if (file.IntArray[i, j, k] == p)
                            {
                                Histo[i, p] = Histo[i, p] + 1;
                                break;
                            }
                        }
                    }

            //计算累计直方图
            int[,] TempAccuHisto = new int[7, 256];
            double[,] AccuHisto = new double[7, 256];
            for (int i = 0; i < file.Head.Bands; i++)
                for (int p = 0; p < 256; p++)
                {
                    TempAccuHisto[i, p] = 0;
                    AccuHisto[i, p] = 0.00;
                }
            for (int i = 0; i < file.Head.Bands; i++)
            {

                for (int p = 0; p < 256; p++)
                {
                    // TempAccuHisto[i,p]+=Histo[i,p];  为什么错了
                    for (int s = 0; s < p; s++)
                    {
                        TempAccuHisto[i, p] = TempAccuHisto[i, p] + Histo[i, s];
                    }

                    //AccuHisto[i, p] = TempAccuHisto[i, p] /360000;
                    //double xxx = TempAccuHisto[i, p] / 360000;
                    //AccuHisto[i, p] = xxx;
                }
            }

            for (int i = 0; i < file.Head.Bands; i++)
                for (int p = 0; p < 256; p++)
                {
                    AccuHisto[i, p] = 1.000 * TempAccuHisto[i, p] / 360000;//真要吐血了，看了三个小时
                }

            //计算百分之二拉伸的左右区间端点值
            int[,] range = new int[7, 2];

            for (int i = 0; i < file.Head.Bands; i++)
            {

                double mini = AccuHisto[i, 0] - 0.02;  //存储离0.02最近的灰度值
                if (mini < 0)
                    mini = -1 * mini;
                double maxi = AccuHisto[i, 0] - 0.98;  //存储离0.98最近的灰度值
                if (maxi < 0)
                    maxi = -1 * maxi;

                for (int p = 0; p < 256; p++)
                {
                    double a = AccuHisto[i, p] - 0.02;
                    if (a < 0)
                        a = -1 * a;
                    double b = AccuHisto[i, p] - 0.98;
                    if (b < 0)
                        b = -1 * b;
                    if (a < mini)
                    {
                        mini = a;
                        range[i, 0] = p;
                    }
                    if (b < maxi)
                    {
                        maxi = b;
                        range[i, 1] = p;
                    }
                }
            }
            
            //将左右端点值带入拉伸公式
            for (int i = 0; i < file.Head.Bands; i++)
            {
                for (int j = 0; j < file.Head.Lines; j++)
                {
                    for (int k = 0; k < file.Head.Samples; k++)
                    {
                        if (file.IntArray[i, j, k] < range[i, 0])
                        {
                            outFile.IntArray[i, j, k] = 0;
                        }

                        else if (file.IntArray[i, j, k] > range[i, 1])
                        {
                            outFile.IntArray[i, j, k] = 255;
                        }
                        else
                        {
                            outFile.IntArray[i, j, k] = (int)(255 * (file.IntArray[i, j, k] - range[i, 0]) / (range[i, 1] - range[i, 0]));
                        }
                    }
                }
            }

            return outFile;
        }

        // 直方图均衡化      
        public Files Equalization(Files file)
        {
            Files outFile = file.Clone();
            outFile.Name += "_Equalized";

            funcStatistics funcSt = new funcStatistics();
            int[,] frestr = funcSt.Frequency(file);

            int h; double q;
            for (int i = 0; i < file.Head.Bands; i++)
            {
                for (int j = 0; j < file.Head.Lines; j++)
                {
                    for (int k = 0; k < file.Head.Samples; k++)
                    {
                        h = 0;
                        for (int w = 1; w <= file.IntArray[i, j, k]; w++)
                        {
                            h = h + frestr[i, w];
                        }
                        q = 255 * h / (file.Head.Lines * file.Head.Samples * 1.00);
                        outFile.IntArray[i, j, k] = (int)q;
                    }
                }
            }
            return outFile;
        }

        // 直方图规定化
        public Files Specification(Files raw, Files refer)
        {
            funcStretch funcSt = new funcStretch();

            //原始图像
            Files old = raw.Clone();

            //参考图像
            Files temp = refer.Clone();

            Files fir = funcSt.Equalization(old);
            Files sec = funcSt.Equalization(temp);

            int[,] h1 = new int[10, 256];
            int[,] h2 = new int[10, 256];
            int[,] first = new int[10, 256];
            int[,] second = new int[10, 256];
            int[,] final = new int[10, 256];
            int[,,] chayi = new int[10, 256, 256];
            int[,,] mark = new int[10, 1000, 1000];

            for (int i = 0; i < temp.Head.Bands; i++)
            {
                for (int m = 0; m < 256; m++)
                {
                    for (int j = 0; j < temp.Head.Lines; j++)
                    {
                        for (int k = 0; k < temp.Head.Samples; k++)
                        {
                            if (fir.IntArray[i, j, k] == m)
                            {
                                h1[i, m]++;
                            }
                            if (sec.IntArray[i, j, k] == m)
                            {
                                h2[i, m]++;
                            }
                        }
                    }
                }//灰度直方图
            }

            int t1 = 0, t2 = 0;
            for (int h = 0; h < 10; h++)
            {
                for (int n = 0; n < 256; n++)
                {
                    t1 = t1 + h1[h, n];
                    first[h, n] = t1;

                    t2 = t2 + h2[h, n];
                    second[h, n] = t2;
                }
                t1 = 0;
                t2 = 0;
            }//累计直方图

            for (int x = 0; x < 10; x++)
            {
                for (int y = 0; y < 256; y++)
                {
                    for (int z = 0; z < 256; z++)
                    {
                        chayi[x, y, z] = Math.Abs(first[x, y] - second[x, z]);
                    }
                }
            }//两者累计直方图的差异

            int[,] jilu = new int[10, 256];
            int minnum, minvalue;
            for (int d = 0; d < 10; d++)
            {
                for (int g = 0; g < 256; g++)
                {
                    minnum = 0;
                    minvalue = chayi[d, g, 0];//假设1的灰度=g和2灰度=0的累计频数差异最小
                    for (int f = 0; f < 256; f++)
                    {
                        if (minvalue > chayi[d, g, f])
                        {
                            minvalue = chayi[d, g, f];//如果1的灰度=g和2的灰度=f的累计频数差异最小
                            minnum = f;//记录波段2的灰度值
                        }
                    }
                    jilu[d, g] = minnum;//g为波段一的灰度，minnum为与之差异最小的波段2的灰度
                }
            }//记录最小差异

            for (int i = 0; i < temp.Head.Bands; i++)
            {
                for (int a4 = 0; a4 < 256; a4++)
                {
                    for (int j = 0; j < temp.Head.Lines; j++)
                    {
                        for (int k = 0; k < temp.Head.Samples; k++)
                        {
                            if (mark[i, j, k] != 1)
                            {
                                if (fir.IntArray[i, j, k] == a4)
                                {
                                    fir.IntArray[i, j, k] = jilu[i, a4];
                                    mark[i, j, k] = 1;
                                }

                            }

                        }

                    }
                }
            }

            for (int q = 0; q < temp.Head.Bands; q++)
            {
                for (int j = 0; j < temp.Head.Lines; j++)
                {
                    for (int k = 0; k < temp.Head.Samples; k++)
                    {
                        temp.IntArray[q, j, k] = fir.IntArray[q, j, k];
                    }
                }
            }

            return temp;
        }

        // OIF
        public double CalOIF(Files file, int a, int b, int c)
        {
            ////计算7个标准差，直接从RSProject复制而来。  
            //计算相关系数矩阵，如上复制而来
            double[,] Cov = new double[file.Head.Bands, file.Head.Bands];
            double[] ex = new double[file.Head.Bands];
            double[] ex2 = new double[file.Head.Bands];
            double[] dx = new double[file.Head.Bands];//这里指标准差
            double[,] Cor = new double[file.Head.Bands, file.Head.Bands];//相关系数

            //初始化
            for (int i = 0; i < file.Head.Bands; i++)
            {
                for (int j = 0; j < file.Head.Bands; j++)
                {
                    Cov[i, j] = 0;
                }
                ex[i] = 0;
                ex2[i] = 0;
            }

            //计算标准差
            for (int i = 0; i < file.Head.Bands; i++)
            {
                for (int j = 0; j < file.Head.Lines; j++)
                    for (int k = 0; k < file.Head.Samples; k++)
                    {
                        ex[i] += file.IntArray[i, j, k];
                        ex2[i] += Math.Pow(file.IntArray[i, j, k], 2);
                    }
                ex[i] /= (file.Head.Lines * file.Head.Samples);
                ex2[i] /= (file.Head.Lines * file.Head.Samples);
                dx[i] = Math.Sqrt(ex2[i] - Math.Pow(ex[i], 2));
            }

            //计算协方差
            for (int m = 0; m < file.Head.Bands; m++)
            {
                for (int n = 0; n < file.Head.Bands; n++)
                {
                    for (int j = 0; j < file.Head.Lines; j++)
                    {
                        for (int k = 0; k < file.Head.Samples; k++)
                        {
                            Cov[m, n] += (file.IntArray[m, j, k] - ex[m]) * (file.IntArray[n, j, k] - ex[n]);
                        }
                    }
                    Cov[m, n] /= (file.Head.Lines * file.Head.Samples);
                }
            }

            //计算相关系数矩阵
            for (int i = 0; i < file.Head.Bands; i++)
                for (int j = 0; j < file.Head.Bands; j++)
                {
                    Cor[i, j] = Cov[i, j] / (dx[i] * dx[j]);
                }

            double midiup = dx[a + 1] + dx[b + 1] + dx[c + 1];
            if (Cor[a, b] < 0)
                Cor[a, b] = -1 * Cor[a, b];
            if (Cor[a, c] < 0)
                Cor[a, c] = -1 * Cor[a, c];
            if (Cor[c, b] < 0)
                Cor[c, b] = -1 * Cor[c, b];
            double mididown = Cor[a, b] + Cor[a, c] + Cor[b, c];
            double midi = 1.00 * midiup / mididown;

            return midi;
        }

        // 重分级
        public Files Reclass(Files file)
        {
            Files outFile = file.Clone();
            for (int i = 0; i < file.Head.Bands; i++)
            {
                double max = file.IntArray[0, 0, 0], min = file.IntArray[0, 0, 0];
                for (int j = 0; j < file.Head.Lines; j++)
                {
                    for (int k = 0; k < file.Head.Samples; k++)
                    {
                        if (file.IntArray[i, j, k] < min)
                            min = file.IntArray[i, j, k];
                        if (file.IntArray[i, j, k] > max)
                            max = file.IntArray[i, j, k];
                    }
                }
                for (int j = 0; j < file.Head.Lines; j++)
                {
                    for (int k = 0; k < file.Head.Samples; k++)
                    {
                        outFile.IntArray[i, j, k] = 255 * (file.IntArray[i, j, k] - min) / (max - min);
                    }
                }
            }

            return outFile;
        }
    }
}

